using System;

using Atomic.Libraries.Mathematics.LinearAlgebra.AlgLib;

namespace Atomic.Libraries.Mathematics.LinearAlgebra
{
	[Serializable]
	public sealed class EigenvalueDecomposition
	{
		private Vector eigenvalues;
		private Matrix eigenvectors;

		private EigenvalueDecomposition()
		{
		}

		/// <summary>
		/// Decomposes a symmetric real matrix. Only the upper triangular part is used; symmetry is not checked.
		/// </summary>
		public static EigenvalueDecomposition Decompose(Matrix a)
		{
			EigenvalueDecomposition ed = new EigenvalueDecomposition();
			Decompose(a, out ed.eigenvalues, out ed.eigenvectors);
			return ed;
		}

		/// <summary>
		/// Decomposes a symmetric real matrix. Only the upper triangular part is used; symmetry is not checked.
		/// </summary>
		public static EigenvalueDecomposition Decompose(double[,] a)
		{
			return Decompose(new Matrix(a));
		}

		/*public static EigenvalueDecomposition DecomposeSymmetrized(Matrix a)
		{
			if (a.Rows != a.Columns)
			{
				throw new ArgumentException("Matrix isn't a square matrix.");
			}

			return Decompose(0.5 * (a + Matrix.Transpose(a)));
		}*/

		/// <summary>
		/// Decomposes a symmetric real matrix. Only the upper triangular part is used; symmetry is not checked. Eigenvalues are sorted in ascending order. Eigenvectors are stored as matrix columns.
		/// </summary>
		public static void Decompose(Matrix a, out Vector eigenvalues, out Matrix eigenvectors)
		{
			int n = a.Rows;

			if (a.Columns != n)
			{
				throw new ArgumentException("Matrix isn't a square matrix.");
			}

			double[,] a0 = new double[n, n];
			for (int i = 0; i < n; i++)
			{
				for (int j = i; j < n; j++)
				{
					a0[i, j] = a[i, j];
				}
			}

			double[] d = new double[n];
			double[,] z = new double[n, n];

			evd.smatrixevd(a0, n, 1, true, ref d, ref z);

			eigenvalues = new Vector(d);
			eigenvectors = new Matrix(z);
		}

		/// <summary>
		/// Eigenvalues are sorted in ascending order.
		/// </summary>
		public Vector Eigenvalues
		{
			get
			{
				return eigenvalues;
			}
		}

		/// <summary>
		/// Eigenvectors are stored as matrix columns.
		/// </summary>
		public Matrix Eigenvectors
		{
			get
			{
				return eigenvectors;
			}
		}
	}
}
